<html>
 <head>
  <meta charset="utf-8"/>
  <meta content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no" name="viewport"/>
  <title>
   主题：多重比较可视化求助  | 数螺 | NAUT IDEA
  </title>
  <link href="http://cdn.bootcss.com/bootstrap/3.3.6/css/bootstrap-theme.min.css" rel="stylesheet"/>
  <link href="http://cdn.bootcss.com/bootstrap/3.3.6/css/bootstrap.min.css" rel="stylesheet"/>
  <style type="text/css">
   #xmain img {
                  max-width: 100%;
                  display: block;
                  margin-top: 10px;
                  margin-bottom: 10px;
                }

                #xmain p {
                    line-height:150%;
                    font-size: 16px;
                    margin-top: 20px;
                }

                #xmain h2 {
                    font-size: 24px;
                }

                #xmain h3 {
                    font-size: 20px;
                }

                #xmain h4 {
                    font-size: 18px;
                }


                .header {
	           background-color: #0099ff;
	           color: #ffffff;
	           margin-bottom: 20px;
	        }

	        .header p {
                  margin: 0px;
                  padding: 10px 0;
                  display: inline-block;  
                  vertical-align: middle;
                  font-size: 16px;
               }

               .header a {
                 color: white;
               }

              .header img {
                 height: 25px;
              }
  </style>
  <script src="http://cdn.bootcss.com/jquery/3.0.0/jquery.min.js">
  </script>
  <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript">
   MathJax.Hub.Config({elements: ["bbpress-forums"]});
  </script>
  <script src="http://nautstatic-10007657.file.myqcloud.com/static/css/readability.min.js" type="text/javascript">
  </script>
  <script type="text/javascript">
   $(document).ready(function() {
                 var loc = document.location;
                 var uri = {
                  spec: "http://cos.name/cn/topic/142002/?new=1",
                  host: "http://cos.name",
                  prePath: "http://cos.name",
                  scheme: "http",
                  pathBase: "http://cos.name/"
                 };
    
                 var documentClone = document.cloneNode(true);
                 var article = new Readability(uri, documentClone).parse();
     
                 document.getElementById("xmain").innerHTML = article.content;
                });
  </script>
  <!-- 1466445530: Accept with keywords: (title(0.0):主题,可视化,论坛, topn(0.533333333333):方差,画图,帖子,分布,会员,计量,直观,结果,均值,参考文献,可视化,方法,正态分布,条形图,置信区间,函数,公式,差异,问题,普通,分位,误差,置信度,数据,主题,算法,数值,显著,作者,文献).-->
 </head>
 <body class="topic bbpress single single-topic postid-142002 single-author sidebar" onload="">
  <div class="header">
   <div class="container">
    <div class="row">
     <div class="col-xs-6 col-sm-6 text-left">
      <a href="/databee">
       <img src="http://nautidea-10007657.cos.myqcloud.com/logo_white.png"/>
      </a>
      <a href="/databee">
       <p>
        数螺
       </p>
      </a>
     </div>
     <div class="hidden-xs col-sm-6 text-right">
      <p>
       致力于数据科学的推广和知识传播
      </p>
     </div>
    </div>
   </div>
  </div>
  <div class="container text-center">
   <h1>
    主题：多重比较可视化求助
   </h1>
  </div>
  <div class="container" id="xmain">
   <div class="hfeed site" id="page">
    <header class="site-header" id="masthead" role="banner">
     <div id="cos-logo">
      <a href="http://cos.name/cn">
       <img src="http://cos.name/cn/wp-content/themes/COS-forest/images/headers/cos-logo.png"/>
      </a>
     </div>
     <div class="navbar" id="navbar">
      <nav class="navigation main-navigation" id="site-navigation" role="navigation">
       <h3 class="menu-toggle">
        菜单
       </h3>
       <div class="menu-%e8%8f%9c%e5%8d%951-container">
        <ul class="nav-menu" id="menu-%e8%8f%9c%e5%8d%951">
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-home menu-item-407772" id="menu-item-407772">
          <a href="http://cos.name/cn/">
           论坛首页
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-407773" id="menu-item-407773">
          <a href="http://cos.name/cn/forums/">
           讨论区
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-407774" id="menu-item-407774">
          <a href="http://cos.name/cn/wp-login.php?action=register">
           注册
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-407819" id="menu-item-407819">
          <a href="http://cos.name/">
           主站
          </a>
         </li>
        </ul>
       </div>
      </nav>
      <!-- #site-navigation -->
     </div>
     <!-- #navbar -->
    </header>
    <!-- #masthead -->
    <div class="site-main" id="main">
     <div class="content-area" id="primary">
      <div class="site-content" id="content" role="main">
       <article class="post-142002 topic type-topic status-publish hentry" id="post-142002">
        <header class="entry-header">
         <h1 class="entry-title">
          多重比较可视化求助
         </h1>
        </header>
        <!-- .entry-header -->
        <div class="entry-content">
         <div id="bbpress-forums">
          <div class="bbp-breadcrumb">
           <p>
            <a class="bbp-breadcrumb-home" href="http://cos.name/cn/">
             COS论坛 | 统计之都
            </a>
            <span class="bbp-breadcrumb-sep">
             ›
            </span>
            <a class="bbp-breadcrumb-root" href="http://cos.name/cn/forums/">
             讨论区
            </a>
            <span class="bbp-breadcrumb-sep">
             ›
            </span>
            <a class="bbp-breadcrumb-forum" href="http://cos.name/cn/forum/software/">
             软件应用
            </a>
            <span class="bbp-breadcrumb-sep">
             ›
            </span>
            <a class="bbp-breadcrumb-forum" href="http://cos.name/cn/forum/software/r-language/">
             S-Plus &amp; R语言
            </a>
            <span class="bbp-breadcrumb-sep">
             ›
            </span>
            <span class="bbp-breadcrumb-current">
             多重比较可视化求助
            </span>
           </p>
          </div>
          <div class="bbp-template-notice info">
           <p class="bbp-topic-description">
            该主题包含 19 条回复，6个帖子，最后由
            <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
             <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=14&amp;d=monsterid&amp;r=g"/>
            </a>
            <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
             yufree
            </a>
            在
            <a href="http://cos.name/cn/topic/142002/page/2/#post-389950" title="回复：多重比较可视化求助">
             2 年, 5 月 之前
            </a>
            更新。
           </p>
          </div>
          <div class="bbp-pagination">
           <div class="bbp-pagination-count">
            查看 15 个帖子 - 1 到 15（总计 20 个）
           </div>
           <div class="bbp-pagination-links">
            <span class="page-numbers current">
             1
            </span>
            <a class="page-numbers" href="http://cos.name/cn/topic/142002/page/2/?new=1">
             2
            </a>
            <a class="next page-numbers" href="http://cos.name/cn/topic/142002/page/2/?new=1">
             →
            </a>
           </div>
          </div>
          <ul class="forums bbp-replies" id="topic-142002-replies">
           <li class="bbp-header">
            <div class="bbp-reply-author">
             作者
            </div>
            <!-- .bbp-reply-author -->
            <div class="bbp-reply-content">
             帖子
            </div>
            <!-- .bbp-reply-content -->
           </li>
           <!-- .bbp-header -->
           <li class="bbp-body">
            <div class="bbp-reply-header" id="post-142002">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月11日 上午11:26
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-142002">
               1 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-142002 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-1 user-id-379087 topic-author post-142002 topic type-topic status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               近期处理一批数据，三组，方差不齐，个数不一, 要求进行方差分析与多重比较并要求把多重比较结果用条形图画出来。
              </p>
              <p>
               分析过程（可忽略，直接看问题）
              </p>
              <p>
               条形图上加个误差线不难，两两比较给出置信区间也不难，但加多重比较结果叠加到反应均值的条形图上好像只能手工添加。查了下发现有一种多重比较的方法Gabriel comparison interval(Gabriel, 1978)比较牛，直接出图就OK。不过这是在SAS里实现的，本着厚脸皮偷懒的心到网上搜R的函数，但没找到，转念想想不就是算出个区间加到条形图上，自己来就是了。然后找到了Gabriel, K.R. 1978. A simple method of multiple comparison of means. J. Amer. Stat. Assoc. 73: 724-729.的原文，读了下就泄气了。别的都好说，但Studentized Maximum Modulus要查表，估计这年代守着电脑查表不够水平，就打算找找这个Studentized Maximum Modulus的分布，然后根据置信度，自由度与水平数（在我这里就是3）得到这个统计量。找了半天只找到一个IBM网站上的公式：
              </p>
              <p>
               <a href="http://pic.dhe.ibm.com/infocenter/spssstat/v20r0m0/index.jsp?topic=%2Fcom.ibm.spss.statistics.help%2Falg_cdf-pdf-rand_continuous_smod.htm" rel="nofollow">
                http://pic.dhe.ibm.com/infocenter/spssstat/v20r0m0/index.jsp?topic=%2Fcom.ibm.spss.statistics.help%2Falg_cdf-pdf-rand_continuous_smod.htm
               </a>
              </p>
              <p>
               然后问题就在这里，这个公式里有个gamma分布函数，后面有个2，我搞不清楚这是乘法还是平方。然后又搜了半天，由于非统计科班，找了半天也没找到跟IBM长得像的函数。之后我打算上stackoverflow碰碰运气，然后就被华丽的鄙视了。鄙视过程参见http://stackoverflow.com/questions/20475092/how-to-plot-gabriel-comparison-interval-in-r
              </p>
              <p>
               后来又翻了下stats包发现有一个tukey分布，我猜测这个tukey分布跟我要找的Studentized Maximum Modulus有关系，希望这一点得到统计专业的指点。
              </p>
              <p>
               再后来发现Hochberg’s GF2方法也用了Studentized Maximum Modulus这个统计量，就打算从已有的R包源函数里把这部分翻出来，然后我又天真了，没找到这个R函数。又看到Dunnett-Tukey-Kramer方法可用在方差与个数不同时的检验，然后翻这个DTK包的源函数，又发现了tukey分布的身影
              </p>
              <p>
              </p>
              <pre class="highlight ">SR = qtukey(p = a, nmeans = k, df = df, lower.tail = FALSE)</pre>
              <p>
              </p>
              <p>
               现在对比原文献如果tukey给出的就是m值，那么剩下的工作就轻松加愉快了，把DTK的函数改一下就得到Gabriel的方法，加工一下就可以画图了，甚至都不用纠结于Gabriel方法，直接用DTK方法的结果去画图。具体就是先按照均值排列条形图，这样相邻的误差线就选取两者差值离0最近的值，然后就是看上限与下限是否重叠，不重叠有差异，重叠无差异。
              </p>
              <p>
               如果将多重比较的结果与均值放到一个图里，基本的思路通过误差线表示差异显著性，有重叠就是不显著，清晰直观，这应该是Gabriel的最初思路，个人感觉也很实用，避免了通过条形图均值标准差标准误或95区间去推断多重比较引发的误会。但有两个问题绕不过去如下：
              </p>
              <p>
               问题
              </p>
              <p>
               1 tukey分布可以用来求Studentized Maximum Modulus吗？
              </p>
              <p>
               2 如果不用Gabriel方法做图，用TukeyHSD（假设个数一样），那么可视化过程的误差线选哪一对，两两差异是否会导致与其他对比的显著差异体现不出（如ABC三组，AB差异显著，BC差异显著，AC差异不显著，是否可能BC的下限体现不出AC的差异，ABC均值1，2，3.AB差异1+-0.5,BC差异1+-0.3，AC差异2+-2.5，这时C的下误差线如何选取来体现差异）或者说理论上不可能出现这种情况。有没有更好的解决思路。
              </p>
              <p>
               注：因为不是统计学专业，很多论文没权限，现在看这些概念如雾（di）里（du）看花，请重喷指出不足。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-384973">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月12日 上午9:23
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-384973">
               2 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-384973 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-2 user-id-379087 topic-author post-384973 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               重读了原始文献，自问自答吧。
              </p>
              <p>
               1 两个不是一回事，这里写的很清楚http://pic.dhe.ibm.com/infocenter/spssstat/v20r0m0/index.jsp?topic=%2Fcom.ibm.spss.statistics.help%2Falg_posthoc_range-maxmodulus.htm
              </p>
              <p>
               2 不用Gabriel的计算方法是不行的，SNK也好，Hochberg’s GT2也好，这些多重比较都涉及两两间比较，无法想像在一张图上用一条误差线描述多个比较的结果。
              </p>
              <p>
               不过Gabriel的方法很有启发性，他用一条误差线表示出了所有对比的差异，只不过这个差异符合SMM分布而已。自然是用了很保守的态度，置信度会超过alpha。这样作出的l-u图可用来直观对比差异但必须要注明这个图误差线是表示对比而不是均值描述，而相应的对均值的描述就牺牲掉了，否则还是让人看不懂。
              </p>
              <p>
               虽然样本数不相同时要用到Studentized Maximum Modulus这个统计量（接近GT2方法），但样本量相同时，该方法可用TukeyHSD的q值等效替代。
              </p>
              <p>
               也就是说，先利用
              </p>
              <p>
              </p>
              <pre class="highlight ">SR = qtukey(p = a, nmeans = k, df = df, lower.tail = FALSE)</pre>
              <p>
              </p>
              <p>
               求出q值，再根据Gabriel论文中的公式就可以做图了。至于样本数不相等，就只能通过求解SMM的积分函数来得到M值，由于R中没内置这个分布，我查到的公式不够清晰，所以也就没办法了（其实是不会求解，而这个问题也被晾了好久了http://r.789695.n4.nabble.com/Studentized-maximum-distribution-td836659.html）
              </p>
              <p>
               13日update 网上找到一个SMM的表，现在把Gabriel这个直观的方法写成rgabriel包放到github上了，补充完毕后考虑推到CRAN上。欢迎大家试用https://github.com/yufree/rgabriel
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-386086">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月16日 上午7:09
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-386086">
               3 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-386086 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-3 user-id-379087 topic-author post-386086 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               该R包已经发布在CRAN上了
               <a href="http://cran.r-project.org/web/packages/rgabriel/index.html" rel="nofollow">
                http://cran.r-project.org/web/packages/rgabriel/index.html
               </a>
              </p>
              <p>
               需要注明的是画出的柱形图只用来作多重比较（且方法限于Gabriel（1978）的文章中的方法，你若认可其可能增大的I型错误率就没什么可说的了），并不表示均值的测量误差什么的，当然可以用那个画图函数去手动添加，设定好upper跟lower就可以了。
              </p>
              <p>
               下一步计划是积分，心中对SMM分布的那个表耿耿于怀，现在只能求置信度0.01跟0.05两个，多种比较也不能超过8组，但实在是对R中数值求解没概念，也没找到已知概率密度函数的参数反求SMM值的方法，先这样吧，有精力再处理。
              </p>
              <p>
               ps 赞rstudio的R包开发功能，上手简单，一条龙服务啊。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-386096">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月16日 上午7:24
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-386096">
               4 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-386096 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-4 user-id-376892 post-386096 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/376892/" rel="nofollow" title="查看masterlxh的档案">
               <img src="http://sdn.geekzu.org/avatar/5234856bdacc8486ceb94c69ddd105f6?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/376892/" rel="nofollow" title="查看masterlxh的档案">
               masterlxh
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               tukey老爷子的那本蓝皮的《探索性数据分析》不知道大家还记得不？
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388131">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月22日 上午6:08
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388131">
               5 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388131 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-5 user-id-1 post-388131 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               <img src="http://sdn.geekzu.org/avatar/1022d8e6ebc94e8f6bca9a86cebe312a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               谢益辉
              </a>
              <br/>
              <div class="bbp-author-role">
               站长
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               既然没有解析解，学术以及帝都雾都大，在家偷懒敲键盘整个模拟解应该可以吧。反正计算便宜，就多算点好了，比起看什么Legendre 16点逼近应该舒服一些。
              </p>
              <p>
              </p>
              <pre class="highlight "># 生成一个随机数
smm_one = function(m) max(abs(x &lt;- rnorm(m))) / sd(x)
# 生成n个随机数
rsmm = function(n, m) replicate(n, smm_one(m))
# 求分位数
qsmm = function(q, n = 10000, m) quantile(rsmm(n, m), q)
# 累积概率密度函数
psmm = function(x, n = 10000, m) mean(rsmm(n, m) &lt;= x)</pre>
              <p>
              </p>
              <p>
               话说你动作真够快的，呼哧一下就是一个R包上CRAN了，这种顺利程度建议出门买彩票去。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388153">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月22日 下午4:14
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388153">
               6 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388153 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-6 user-id-379087 topic-author post-388153 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               回复 第5楼 的 谢益辉：
              </p>
              <p>
               谢谢大指点，模拟的方法太妙了。
              </p>
              <p>
               我查了下这个分布应该有分位数、两两对比次数跟自由度三个参数，模拟中m应该指代的是两两对比次数，那自由度这个参数好像没体现出来。sd(x)给出的是单次模拟所得标准差，这个方差乘以自由度df除以生成m的正态分布的方差符合一个自由度为df的卡方分布。不考虑自由度得到的值与自由度趋向无穷大的SMM表数据相近。那么如果考虑自由度能否写成下面这样,或者说我想多了：
              </p>
              <p>
              </p>
              <pre class="highlight ">smm_one = function(m,df) max(abs(rnorm(m)-mean(rnorm(m))))/sqrt(rchisq(1,df)/df)
</pre>
              <p>
              </p>
              <p>
               其实这个包不大顺利，因为没经验，第一次提交的包都过不了check，被Kurt Hornik给打回来了，不过大体上比着DTK这个包写的，第二次提交很快就被Brian Ripley放上去了，几个小时后Uwe Ligges就把windows版给提交了，十分惊叹于r-core的效率，前后不过两天。
              </p>
              <p>
               其实最近翻了turkey分布的源文件，发现是c写的，而且也是先查表后计算，作者好像就是Kurt Hornik，因为c早忘干净了，就打消写SMM分布的心思了。不过现在看看也是死心眼，连模拟这个计算机的强项都忘了。
              </p>
              <p>
               ======
              </p>
              <p>
               23日更新
              </p>
              <p>
               昨晚上那个写法有问题，跟表数据对不上，我想来想去就顾名思义了，最大学生化模数考虑自由度，要不就用t分布去模拟了，写成下面这样，跟表格数据对比了下，数值偏小，但随自由度变化了。现在有点凌乱了，究竟表格数据准还是模拟准，没解析解我倾向于模拟更准。
              </p>
              <p>
              </p>
              <pre class="highlight ">smm_one = function(m,df) max(abs(x &lt;- rt(m,df)))
</pre>
              <p>
              </p>
              <p>
               =====
              </p>
              <p>
               再更
              </p>
              <p>
               用上面得到数值对比已有的表数据偏大而不是偏小，用谢大的公式偏小，考虑到多重比较容易假阳性，暂时采用上面公式，损失些功效，这样先在github上更新，暂时不推cran，等搞到原始文献连同GT2方法一并更新
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388677">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月25日 上午7:37
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388677">
               7 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388677 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-7 user-id-1 post-388677 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               <img src="http://sdn.geekzu.org/avatar/1022d8e6ebc94e8f6bca9a86cebe312a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               谢益辉
              </a>
              <br/>
              <div class="bbp-author-role">
               站长
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               你后面更新的是对的，我刚开始没有仔细看公式。我重新看了一下相关的文献，这个SMM确实应该是t分布的绝对值，随机数可以这样生成：
              </p>
              <p>
              </p>
              <pre class="highlight ">rsmm = function(n, m, df) apply(matrix(abs(rt(m * n, df)), m, n), 2, max)</pre>
              <p>
              </p>
              <p>
               不小心被你带到这条沟里来了，索性看了一下这篇文章：http://www.jstor.org/discover/10.2307/1268584 用他给的公式(2.2)暴力算了一下他的表1，结果能吻合：
              </p>
              <p>
              </p>
              <pre class="highlight ">psmm_x = function(x, c, r, nu) {
  snu = sqrt(nu)
  sx  = snu * x  # for the scaled Chi distribution
  lgx = log(snu) - lgamma(nu/2) + (1 - nu/2) * log(2) + (nu - 1) * log(sx) + (-sx^2/2)
  exp(r * log(2 * pnorm(c * x) - 1) + lgx)
}
psmm = function(x, r, nu) {
  res = integrate(psmm_x, 0, Inf, c = x, r = r, nu = nu)
  res$value
}
qsmm = function(q, r, nu) {
  r = (r * (r - 1)/2)
  if (!is.finite(nu)) return(qnorm(1 - .5 *(1 - q^(1/r))))
  res = uniroot(function(c, r, nu, q) {
    psmm(c, r = r, nu = nu) - q
  }, c(0, 100), r = r, nu = nu, q = q)
  res$root
}
qsmm(.99, r=3, nu=5)
qsmm(.99, r=3, nu=Inf)
ks = 3:20
vs = c(5, 7, 10, 12, 16, 20, 24, 30, 40, 60, 120, Inf)
q0.01 = matrix(nrow = length(ks), ncol = length(vs), dimnames = list(ks, vs))
for (k in ks) {
  for (v in vs) {
    q0.01[as.character(k), as.character(v)] = qsmm(.99, r = k, nu = v)
  }
}
round(q0.01, 3)</pre>
              <p>
              </p>
              <p>
               代码写得比较混乱，大概意思是用uniroot()求F(x) = q的根（即分位数），累积概率分布函数是用积分integrate()求的，被积函数为了数值稳定性先取log()再exp()。当自由度无穷大时，t分布就是正态分布，文中的公式(3.5)给出了分位数算法，就不用算积分那么麻烦了。
              </p>
              <p>
               虽然结果能吻合，但还是疑云重重，因为毕竟模拟出来的结果和他给的结果还是有差异的，有时候大，有时候小。说实话我有点怀疑那个公式有问题，不过作者说是从参考文献5和8里来的，而我找到了5，发现里面没有他说的公式2.2（童话里都是骗银滴），文献8是一本1966年的书，书这种东西比较难找，我挖坑只能挖到这里了，接下来把锹还给你接着挖吧，搞不好能发篇文章正本清源发现大家用了四十多年的多重比较是错的呢。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388736">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月25日 下午2:27
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388736">
               8 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388736 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-8 user-id-368512 post-388736 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/368512/" rel="nofollow" title="查看doctorjxd的档案">
               <img src="http://sdn.geekzu.org/avatar/eb9a4b908da88cc2ca797065f49c0f4b?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/368512/" rel="nofollow" title="查看doctorjxd的档案">
               doctorjxd
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               牛！
               <br/>
               就一个字。
               <br/>
               楼主和谢大都是学术极客啊。
               <br/>
               期待早日推翻那个万恶的多重比较方法。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388738">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月25日 下午4:06
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388738">
               9 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388738 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-9 user-id-379087 topic-author post-388738 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               回复 第7楼 的 谢益辉：
              </p>
              <p>
               再次感谢谢大指点！
              </p>
              <p>
               看来这次真的挖了大坑了，目前正在求文献，上面提到的文献我都下不下来准备周末去国科图查查。1966年的书我猜测是Rupert G. Miller, Jr的《Simultaneous Statistical Inference》的前一版，这本书81年第二版有个章节Erratum to: Addendum New Table of the Studentized Maximum Modulus，看意思曾经进行过勘误，但从序言中看到用的是71年的数据，因为还是看不到原文，等我周末把这个疑点去掉。
              </p>
              <p>
               但如果79年用的公式与71年与66年一致却又跟模拟结果不一致，那就真的有意思了，这说明最初的公式可能有问题。可以设想使用最大的t值分布的初衷更多是从直观理解角度，有了想法然后写出其数学描述。后人学这个方法更多关心是直观表述，理解了其设计统计量的初衷剩下的交给计算机或者工具表，如果最初分布描述就错了可能真就长时间没人意识到。我看到sas介绍smm分布也引用了这几篇文献，spss可能也是，这样我更倾向于在R中直接去掉表全用模拟的方法求这个值了，因为这样最符合Gabriel原始文献的本意，基于smm分布的还有GT2跟Dunnett’s T3这两种多重比较的方法。甚至都可以验证下用Studentized Range的HSD方法，如果这两个统计量分布的数学描述都是一个思路写出来的，那很有可能都有问题（更新：看公式应该不是一个路子出来的）。
              </p>
              <p>
               同样如果计算与公式没问题，那么模拟出现的不稳定就可能需要二次模拟来解决。事实上，在周一验证模拟结果与表的数值时，虽然是生成了10000个数值来模拟分布，但得到的smm值分位数忽大忽小。对此我对同一个分位数模拟了100次生成了柱形图，可能次数比较少（笔记本是6年前的老机子），得到的数取均值会稳定一些，但再把这个取均值重复几次波动还是不小，也发现了模拟跟表格数据忽大忽小的趋势，随机模拟本身例如模拟次数可能也有精度限制，并且可能是分布参数的函数。
              </p>
              <p>
              </p>
              <pre class="highlight ">a &lt;- c(0)
for(i in 1:100){
a[i] &lt;- qsmm(0.975,10000,28,200)
}
mean(a)
</pre>
              <p>
              </p>
              <p>
               年底全是总结实验什么的，忙过这阵我好好想想这些问题总结下这个包引发的坑，不过明显感到自己外行了，力不从心。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388739">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月25日 下午4:07
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388739">
               10 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388739 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-10 user-id-379087 topic-author post-388739 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               额 发重复了 留着周末更新
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388824">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月26日 上午2:33
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388824">
               11 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388824 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-11 user-id-111004 post-388824 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/111004/" rel="nofollow" title="查看统计天下的档案">
               <img src="http://sdn.geekzu.org/avatar/a22f91bab0c7eeb12e67c53af855186a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/111004/" rel="nofollow" title="查看统计天下的档案">
               统计天下
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               怎么上传附件啊，我把这个包做到SPSS中去了。算了，放微博吧：http://weibo.com/1261253730/Ap3cqagSf?ref=home
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-388985">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月26日 上午7:11
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-388985">
               12 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-388985 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-12 user-id-379087 topic-author post-388985 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               回复 第11楼 的 统计天下：
              </p>
              <p>
               谢谢！
              </p>
              <p>
               不过spss在one-way anova里提供了gabriel多重比较，但没提供作图。目前这个包没有按文献的密度分布函数去写，只是提供了表插值，而spss应该是按公式算的。如果那个计算公式没问题，实际上还是调用内置函数里差值的95置信区间去作图会准一些。另外sas里也提供了这个函数，目前坑在那个计算公式了，我是打算下一版中去掉表全用模拟的方法或者先本地模拟生成一个表，不过得等查到原始文献，别折腾半天发现模拟本身会错意了就麻烦了。现在有谢大的暴力算法，完全用公式也不是问题，只是得确认无误才行（如果这样说到底整个算法都是谢大写的，我就是加了个壳，下次更新时会把谢大列为作者而不仅仅是感谢了）。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-389014">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月26日 上午8:14
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-389014">
               13 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-389014 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-13 user-id-111004 post-389014 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/111004/" rel="nofollow" title="查看统计天下的档案">
               <img src="http://sdn.geekzu.org/avatar/a22f91bab0c7eeb12e67c53af855186a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/111004/" rel="nofollow" title="查看统计天下的档案">
               统计天下
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               回复 第12楼 的 yufree：好啊，谢大更新，你更新，我更新。慢着，我好像懂了什么，大概是幼儿园的苹果歌，你一个我一个，好欢乐。[s:11]
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-389514">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月27日 上午4:16
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-389514">
               14 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-389514 -->
            <div class="even bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-14 user-id-1 post-389514 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               <img src="http://sdn.geekzu.org/avatar/1022d8e6ebc94e8f6bca9a86cebe312a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/1/" rel="nofollow" title="查看谢益辉的档案">
               谢益辉
              </a>
              <br/>
              <div class="bbp-author-role">
               站长
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               我提到的两篇文献在这里：
               <br/>
               <a href="https://www.dropbox.com/s/1g9xjzv87t7otdh/1268584.pdf" rel="nofollow">
                https://www.dropbox.com/s/1g9xjzv87t7otdh/1268584.pdf
               </a>
               <br/>
               <a href="https://www.dropbox.com/s/fov1m4mm3tlqv4a/2334520.pdf" rel="nofollow">
                https://www.dropbox.com/s/fov1m4mm3tlqv4a/2334520.pdf
               </a>
              </p>
              <p>
               那本书确实是Simultaneous Statistical Inference。亚马逊上有预览，不过预览里关键的几页看不到，我们图书馆有这本书，我懒得爬出去找了。
              </p>
              <p>
               今天我再想了一下这个问题，意识到它其实比我原先想象的要简单得多，因为我比较懒，能模拟的情况下先考虑模拟，但模拟其实是最劣等的方案，最理想的方案当然是找到解析解，对这个问题来说解析解不太可能。折衷一下，半计算半解析：
              </p>
              <p>
               若X服从t分布，累积分布函数是F(x)，密度函数是f(x)，令Y = |X|，那么密度g(y) = 2f(y)，y &gt; 0；分布G(y) = 2F(y) – 1。对r个独立的t分布随机变量X1, …, Xn，令Z = max(|X1|, …, |Xn|)，那么Z分布H(z) = (2F(z) – 1)^r。这都是些比较初等的概率推导。
              </p>
              <p>
               既然分布函数能以很简单的形式写出来，代码也就大大简化了，可以直接用t分布的分布函数
               <code>
                pt()
               </code>
               ：
              </p>
              <p>
              </p>
              <pre class="highlight ">psmm = function(x, r, nu) {
  ifelse(x &gt; 0, (2 * pt(x, nu) - 1)^r, 0)
}
# 分位数函数仍然用H(x) - q = 0求解x
qsmm = function(q, r, nu) {
  res = uniroot(function(x, r, nu, q) {
    psmm(x, r = r, nu = nu) - q
  }, c(0, 100), r = r, nu = nu, q = q)
  res$root
}
# 算0.99分位数
qsmm(.99, r=3, nu=5)
ks = 3:20
vs = c(5, 7, 10, 12, 16, 20, 24, 30, 40, 60, 120, Inf)
q0.01 = matrix(nrow = length(ks), ncol = length(vs), dimnames = list(ks, vs))
for (k in ks) {
  for (v in vs) {
    q0.01[as.character(k), as.character(v)] = qsmm(.99, r = k*(k-1)/2, nu = v)
  }
}
cbind('r*' = choose(ks, 2), round(q0.01, 3))
    r*      5     7    10    12    16    20    24    30    40    60   120   Inf
3    3  5.243 4.353 3.825 3.647 3.443 3.329 3.257 3.188 3.121 3.056 2.994 2.934
4    6  6.133 4.941 4.256 4.029 3.771 3.629 3.539 3.453 3.370 3.291 3.215 3.143
5   10  6.862 5.404 4.584 4.315 4.013 3.848 3.744 3.644 3.549 3.459 3.372 3.289
6   15  7.491 5.791 4.852 4.547 4.206 4.021 3.905 3.794 3.689 3.589 3.493 3.402
7   21  8.051 6.126 5.079 4.742 4.368 4.165 4.038 3.918 3.803 3.695 3.591 3.493
8   28  8.558 6.424 5.278 4.912 4.506 4.288 4.152 4.023 3.900 3.784 3.673 3.569
9   36  9.024 6.693 5.454 5.061 4.628 4.396 4.251 4.114 3.984 3.861 3.744 3.634
10  45  9.457 6.939 5.614 5.196 4.737 4.491 4.339 4.194 4.058 3.929 3.807 3.691
11  55  9.862 7.166 5.760 5.319 4.835 4.577 4.417 4.267 4.124 3.989 3.862 3.742
12  66 10.244 7.378 5.894 5.431 4.925 4.656 4.489 4.332 4.184 4.044 3.912 3.787
13  78 10.606 7.576 6.019 5.535 5.008 4.728 4.555 4.392 4.238 4.094 3.957 3.829
14  91 10.950 7.762 6.135 5.632 5.084 4.794 4.615 4.447 4.288 4.139 3.999 3.866
15 105 11.279 7.939 6.245 5.722 5.156 4.856 4.672 4.498 4.335 4.181 4.037 3.901
16 120 11.595 8.107 6.348 5.807 5.223 4.914 4.724 4.546 4.378 4.221 4.073 3.933
17 136 11.898 8.267 6.446 5.888 5.286 4.968 4.773 4.590 4.418 4.257 4.106 3.963
18 153 12.190 8.419 6.538 5.964 5.345 5.020 4.820 4.632 4.456 4.291 4.137 3.991
19 171 12.472 8.566 6.627 6.037 5.402 5.068 4.863 4.672 4.492 4.324 4.166 4.018
20 190 12.745 8.707 6.711 6.106 5.455 5.114 4.905 4.709 4.526 4.354 4.193 4.042</pre>
              <p>
              </p>
              <p>
               文献中的分位数表是下面这个，可以看出自由度无穷大的时候能对得上（正态分布），但自由度比较小的时候他的分位数值总体偏小（相应置信区间的长度也会小）：
              </p>
              <p>
               <img src="http://i.imgur.com/uxzspO6.png"/>
              </p>
              <p>
               要么是 a) 参考文献13（Stoline等人）的Fortran程序写的有问题，要么是 b) R的t分布函数写的有问题，要么是 c) 那1966年的书推导有问题。现在看来 a) 和 c) 应该是一致的，因为他们的结果我用暴力积分的办法可以重复出来。那么要么R错了，要么公式错了。希望不是R搞错了，不然感觉真的不会再爱了。
              </p>
              <p>
               20分钟后更新
               <br/>
               ===
              </p>
              <p>
               再看了一遍公式2.2，发现它真的错了耶。
              </p>
              <p>
               <img src="http://i.imgur.com/JNE3Wcz.png"/>
              </p>
              <p>
               那个幂r应该提到积分的外面去。我了个去，我这榆木脑袋居然能发现一处错误的公式，生平第一次哎。叩谢楼主的这份圣诞礼物。
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
            <div class="bbp-reply-header" id="post-389526">
             <div class="bbp-meta">
              <span class="bbp-reply-post-date">
               2013年12月27日 上午6:09
              </span>
              <a class="bbp-reply-permalink" href="http://cos.name/cn/topic/142002/#post-389526">
               15 楼
              </a>
              <span class="bbp-admin-links">
              </span>
             </div>
             <!-- .bbp-meta -->
            </div>
            <!-- #post-389526 -->
            <div class="odd bbp-parent-forum-999 bbp-parent-topic-142002 bbp-reply-position-15 user-id-379087 topic-author post-389526 reply type-reply status-publish hentry">
             <div class="bbp-reply-author">
              <a class="bbp-author-avatar" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               <img src="http://sdn.geekzu.org/avatar/1bfcfcb9cd4c45540ff1462b41c1a01a?s=80&amp;d=monsterid&amp;r=g"/>
              </a>
              <br/>
              <a class="bbp-author-name" href="http://cos.name/cn/profile/379087/" rel="nofollow" title="查看yufree的档案">
               yufree
              </a>
              <br/>
              <div class="bbp-author-role">
               普通会员
              </div>
             </div>
             <!-- .bbp-reply-author -->
             <div class="bbp-reply-content">
              <p>
               回复 第14楼 的 谢益辉：
              </p>
              <p>
               恭喜谢大，我已在github添加谢大为rgabriel包作者，gabriel函数也已经重写，今天要有些着急的杂事，文档与版本的更新稍后处理。
              </p>
              <p>
               看来还应该是内行办内行的事，我看到那个公式也很难想到是哪里错了，期待谢大写paper作为rgabriel包的参考文献。
              </p>
              <p>
               此外，鉴于SAS与SPSS可能使用了Stoline等人的那个算法或者那个表，目前不建议使用这两个软件进行包括Gabriel在内的使用SMM分布的多重比较并建议SAS用户与SPSS用户按照谢大给出的表验证SMM分布（在下穷学生一枚，买不起这两个软件）的正确性。
              </p>
              <p>
               ps：现在对模拟有兴趣了，有时间比较下多次小数模拟与单次大数模拟的对比，看同样精度哪种省CPU。
              </p>
              <p>
               =====
              </p>
              <p>
               大概更新了一下文档，目前0.6版，写的比较凌乱，后续整理好推CRAN。
              </p>
              <p>
              </p>
              <pre class="highlight ">library('devtools')
install_github('rgabriel', 'yufree')
</pre>
              <p>
              </p>
              <p>
               或者
              </p>
              <p>
               <a href="https://dl.dropboxusercontent.com/u/23219669/rgabriel-achieve/rgabriel_0.6.tar.gz" rel="nofollow">
                https://dl.dropboxusercontent.com/u/23219669/rgabriel-achieve/rgabriel_0.6.tar.gz
               </a>
              </p>
             </div>
             <!-- .bbp-reply-content -->
            </div>
            <!-- .reply -->
           </li>
           <!-- .bbp-body -->
           <li class="bbp-footer">
            <div class="bbp-reply-author">
             作者
            </div>
            <div class="bbp-reply-content">
             帖子
            </div>
            <!-- .bbp-reply-content -->
           </li>
           <!-- .bbp-footer -->
          </ul>
          <!-- #topic-142002-replies -->
          <div class="bbp-pagination">
           <div class="bbp-pagination-count">
            查看 15 个帖子 - 1 到 15（总计 20 个）
           </div>
           <div class="bbp-pagination-links">
            <span class="page-numbers current">
             1
            </span>
            <a class="page-numbers" href="http://cos.name/cn/topic/142002/page/2/?new=1">
             2
            </a>
            <a class="next page-numbers" href="http://cos.name/cn/topic/142002/page/2/?new=1">
             →
            </a>
           </div>
          </div>
          <div class="bbp-no-reply" id="no-reply-142002">
           <div class="bbp-template-notice">
            <p>
             您必须先登录才能回复该主题。
            </p>
           </div>
          </div>
         </div>
        </div>
        <!-- .entry-content -->
        <footer class="entry-meta">
        </footer>
        <!-- .entry-meta -->
       </article>
       <!-- #post -->
       <div class="comments-area" id="comments">
       </div>
       <!-- #comments -->
      </div>
      <!-- #content -->
     </div>
     <!-- #primary -->
     <div class="sidebar-container" id="tertiary" role="complementary">
      <div class="sidebar-inner">
       <div class="widget-area">
        <aside class="widget bbp_widget_login" id="bbp_login_widget-2">
         <h3 class="widget-title">
          登录
         </h3>
         <form action="http://cos.name/cn/wp-login.php" class="bbp-login-form" method="post">
          <fieldset>
           <legend>
            登录
           </legend>
           <div class="bbp-username">
            <label for="user_login">
             用户名:
            </label>
           </div>
           <div class="bbp-password">
            <label for="user_pass">
             密码:
            </label>
           </div>
           <div class="bbp-remember-me">
            <label for="rememberme">
             记住用户名
            </label>
           </div>
           <div class="bbp-submit-wrapper">
            <button class="button submit user-submit" id="user-submit" name="user-submit" tabindex="104" type="submit">
             登录
            </button>
           </div>
           <div class="bbp-login-links">
            <a class="bbp-register-link" href="http://cos.name/cn/wp-login.php?action=register" title="注册">
             注册
            </a>
            <a class="bbp-lostpass-link" href="http://cos.name/cn/wp-login.php?action=lostpassword" title="忘记密码">
             忘记密码
            </a>
           </div>
          </fieldset>
         </form>
        </aside>
        <aside class="widget widget_text" id="text-7">
         <h3 class="widget-title">
          搜索
         </h3>
         <div class="textwidget">
          <form action="http://www.google.com/search" id="bbp-search-form" method="get" onsubmit="Gsitesearch(this)" role="search">
           <div>
           </div>
          </form>
          <form id="bbp-search-form-baidu" onsubmit="g(this)" role="search">
           <div>
           </div>
          </form>
         </div>
        </aside>
        <aside class="widget widget_text" id="text-2">
         <h3 class="widget-title">
          新鲜事
         </h3>
         <div class="textwidget">
          <ul>
           <li>
            <a href="http://cos.name/cn/topics/">
             最新帖子
            </a>
           </li>
           <li>
            <a href="http://cos.name/cn/view/popular/">
             最热门主题
            </a>
           </li>
           <li>
            <a href="http://cos.name/cn/view/no-replies/">
             消灭零回复
            </a>
           </li>
          </ul>
         </div>
        </aside>
        <aside class="widget widget_text" id="text-3">
         <h3 class="widget-title">
          RSS订阅
         </h3>
         <div class="textwidget">
          <ul>
           <li>
            <img src="http://cos.name/wp-includes/images/rss.png"/>
            <a href="http://cos.name/cn/topics/feed/">
             所有主题
            </a>
           </li>
           <li>
            <img src="http://cos.name/wp-includes/images/rss.png"/>
            <a href="http://cos.name/cn/forums/feed/">
             所有帖子
            </a>
           </li>
          </ul>
         </div>
        </aside>
       </div>
       <!-- .widget-area -->
      </div>
      <!-- .sidebar-inner -->
     </div>
     <!-- #tertiary -->
    </div>
    <!-- #main -->
    <footer class="site-footer" id="colophon" role="contentinfo">
     <div class="site-info">
      版权所有 © 2014 统计之都 | 由
      <a href="http://wordpress.org/">
       WordPress
      </a>
      构建 | 主题修改自
      <a href="http://wordpress.org/themes/twentythirteen">
       Twenty Thirteen
      </a>
     </div>
     <!-- .site-info -->
    </footer>
    <!-- #colophon -->
   </div>
   <!-- #page -->
  </div>
 </body>
</html>